#=
计算导数
对e^{V}中的项目求导
=#

#using ReverseDiff: jacobian!, JacobianTape
#using Zygote
using ForwardDiff: jacobian, derivative


function gbar_scratch(ss::ScrollSVD{T}) where T
    siz = size(ss.B[end])
    ptr = findfirst(ss.L)
    if isnothing(ptr)
        VL = Diagonal(ones(siz[1]))
        DL = VL
        UL = VL
        UR, DR, VR = ss.F[end].U, Diagonal(ss.F[end].S), ss.F[end].Vt
        #VL, DL, UL = ss.F[end].U, Diagonal(ss.F[end].S), ss.F[end].Vt
        #UR = Diagonal(ones(siz[1]))
        #DR = UR
        #VR = UR
    elseif ptr == 1
        VL, DL, UL = ss.F[ptr].U, Diagonal(ss.F[ptr].S), ss.F[ptr].Vt
        UR = Diagonal(ones(siz[1]))
        DR = UR
        VR = UR
        #VL = Diagonal(ones(siz[1]))
        #DL = VL
        #UL = VL
        #UR, DR, VR = ss.F[1].U, Diagonal(ss.F[1].S), ss.F[1].Vt
    else
        VL, DL, UL = ss.F[ptr].U, Diagonal(ss.F[ptr].S), ss.F[ptr].Vt
        UR, DR, VR = ss.F[ptr-1].U, Diagonal(ss.F[ptr-1].S), ss.F[ptr-1].Vt
    end
    #gtt = inv(Diagonal(ones(siz[1]))+UR*DR*VR*VL*DL*UL)
    #M = inv(UL*UR) + DR*(VR*VL)*DL
    #Fm = svd(M)
    #gtt = inv(Fm.Vt*UL)*inv(Diagonal(Fm.S))*inv(UR*Fm.U)
    #
    #println(diag(DL))
    #println(diag(DR))
    #M = UL*UR + inv(DL)*adjoint(VR*VL)*inv(DR)
    DLS = Diagonal(ones(Float64, siz[1]))
    DLB = Diagonal(ones(Float64, siz[1]))
    DRS = Diagonal(ones(Float64, siz[1]))
    DRB = Diagonal(ones(Float64, siz[1]))
    for i in Base.OneTo(siz[1])
        if DL[i, i] > 1.0
            DLB[i, i] = DL[i, i]
            DLS[i, i] = 1.0
        else
            DLS[i, i] = DL[i, i]
            DLB[i, i] = 1.0
        end
        if DR[i, i] > 1.0
            DRB[i, i] = DR[i, i]
            DRS[i, i] = 1.0
        else
            DRS[i, i] = DR[i, i]
            DRB[i, i] = 1.0
        end
    end
    #
    M = DLS*UL*UR*DRS + inv(DLB)*adjoint(VR*VL)*inv(DRB)
    #M = DL*UL*UR + adjoint(VR*VL)*inv(DR)
    Fm = svd(M, alg=LinearAlgebra.QRIteration())
    ML = UR*DRS*adjoint(Fm.Vt)
    MR = adjoint(Fm.U)*DLS*UL
    btt = ML*inv(Diagonal(Fm.S))*MR
    #println(Fm.S)
    #btt = UR*adjoint(Fm.Vt)*inv(Diagonal(Fm.S))*adjoint(Fm.U)*UL
    #btt = adjoint(Fm.Vt*adjoint(UR))*inv(Diagonal(Fm.S))*adjoint(adjoint(UL)*Fm.U)
    #btt = UR*DR*adjoint(Fm.Vt)*inv(Diagonal(Fm.S))*adjoint(adjoint(UL)*Fm.U)
    #btt = UR*DR*adjoint(Fm.Vt)*inv(Diagonal(Fm.S))*adjoint(Fm.U)*DL*UL
    #M = DL*UL*UR + adjoint(VR*VL)*inv(DR)
    #Fm = svd(M, alg=LinearAlgebra.QRIteration())
    ##println(Fm.S)
    ##btt = UR*adjoint(Fm.Vt)*inv(Diagonal(Fm.S))*adjoint(Fm.U)*UL
    ##btt = adjoint(Fm.Vt*adjoint(UR))*inv(Diagonal(Fm.S))*adjoint(adjoint(UL)*Fm.U)
    #btt = adjoint(Fm.Vt*adjoint(UR))*inv(Diagonal(Fm.S))*adjoint(Fm.U)*DL*UL
    return btt
end


"""
向下传递一个
"""
function down_prog_gbar(bbar::Matrix{T}, bmat::Matrix{T}) where T
    return inv(bmat)*bbar*bmat
end


"""
从 \\bar{B} = - [ \\bar{G} B^-1 ]^T
"""
function bbar_from_gbar(gbar::Matrix{T}, bmat::Matrix{T}) where T
    return -transpose(gbar * inv(bmat))
end


"""
计算B矩阵的导数，包含数值稳定的办法
"""
function bbar_calc(ss::ScrollSVD{T}, allbmats1::Vector{Matrix{T}}) where T
    Nt = length(allbmats1)
    sslen = length(ss.L)
    Ng = Int(Nt//sslen)
    xsiz = size(allbmats1[end])
    Nx = Int(xsiz[1]//2)
    # 求 \bar  B_tau_xy
    #bbar = Array{ComplexF64, 3}(undef, Nt, 2*Nx, 2*Nx)
    #for bi in Base.OneTo(Nt)
    #    gs1 = rawgbar(allbmats1, bi-1)
    #    bba = bbar_from_gbar(gs1, allbmats1[bi])
    #    bbar[bi, :, :] = bba
    #end
    #在观测的阶段所有，ss中所有内容都应在R中
    if any(ss.L)
        throw(error("has L"))
    end
    #计算gs的顺序是1，N, N-1, , N-Ng+2, 然后推一个R进入L
    bbar2 = Array{ComplexF64, 3}(undef, Nt, 2*Nx, 2*Nx)
    for grpi in Base.OneTo(sslen)
        gss = gbar_scratch(ss)
        for ni in Base.OneTo(Ng)
            gi = Nt - Ng*(grpi-1) - ni + 1
            bi = gi+1 > Nt ? 1 : gi+1
            bbar2[bi, :, :] = bbar_from_gbar(gss, allbmats1[bi])
            #println(grpi, " ", gi)
            gss = down_prog_gbar(gss, allbmats1[gi])
        end
        scrollR2L(ss, [ss.B[sslen-grpi+1]])
        #println(sslen-grpi+1, " ", ss.L)
    end
    #for bi in Base.OneTo(Nt)
    #    println(sum(abs.(bbar[bi, :, :]-bbar2[bi, :, :])))
    #end
    #然后将所有的内容返回R
    for grpi in Base.OneTo(sslen)
        scrollL2R(ss, [ss.B[grpi]])
        #println(grpi, " ", ss.L)
    end
    return bbar2
end

#=
"""
计算B矩阵的正向导数
"""
function pbpn_calc_(hk, hscfg::Matrix{Int}, U::Float64, nx::Vector{Float64}, ny::Vector{Float64}, nz::Vector{Float64}; Z2=true)
    Nt, Nx = size(hscfg)
    #hscfg_f64 = Float64.(hscfg)
    #ReverseDiff
    #fbr(cfg, ny, nz) = bmat_IsingAD(ehk, cfg, U*ones(Nx), zeros(Nx), ny, nz; Z2=Z2)[1]
    #fbi(cfg, ny, nz) = bmat_IsingAD(ehk, cfg, U*ones(Nx), zeros(Nx), ny, nz; Z2=Z2)[2]
    #jfbr = JacobianTape(fbr, (rand(Float64, Nx), rand(Float64, Nx), rand(Float64, Nx)))
    #jfbi = JacobianTape(fbi, (rand(Float64, Nx), rand(Float64, Nx), rand(Float64, Nx)))
    #undef会有奇怪的问题，以后排查一下
    #这里内存用的太多了，优化成随bi和xi走的，节约内存
    #pbrpny = zeros(Nt, 2*Nx, 2*Nx, Nx)#Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    #pbipny = zeros(Nt, 2*Nx, 2*Nx, Nx)#Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    #pbrpnz = zeros(Nt, 2*Nx, 2*Nx, Nx)#Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    #pbipnz = zeros(Nt, 2*Nx, 2*Nx, Nx)#Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    #
    dt = 0.1
    ehk_ = exp(-dt*hk)
    if Z2
        ehk2 = kron([1 0; 0 1], ehk_)
    else
        ehk2 = ehk_
    end
    #fbrxp1(ehk, ny, nz) = bmat_IsingADX(ehk, Val(1), U, 0.0, ny, nz)[1]
    #fbrxm1(ehk, ny, nz) = bmat_IsingADX(ehk, Val(-1), U, 0.0, ny, nz)[1]
    #
    #@time begin
    #for bi in Base.OneTo(Nt)
    #    #println("pbpn ", bi)
    #    #result = (zeros(Float64, 2*Nx, 2*Nx, Nx), zeros(Float64, 2*Nx, 2*Nx, Nx), zeros(Float64, 2*Nx, 2*Nx, Nx))
    #    #jacobian!(result, jfbr, (hscfg_f64[bi, :], ny, nz))
    #    #_, brny, brnz = result
    #    #@time begin
    #    #_, brny, brnz = Zygote.jacobian(fbr, hscfg_f64[bi, :], ny, nz)
    #    #_, biny, binz = Zygote.jacobian(fbi, hscfg_f64[bi, :], ny, nz)
    #    #end
    #    #b1 = reshape(brny, 2*Nx, 2*Nx, Nx)
    #    #b2 = reshape(brnz, 2*Nx, 2*Nx, Nx)
    #    #b3 = reshape(biny, 2*Nx, 2*Nx, Nx)
    #    #b4 = reshape(binz, 2*Nx, 2*Nx, Nx)
    #    #println(b1[1, 1:3, 1], b1[2, 1:3, 1], b1[19, 1:3, 1])
    #    #@time begin
    #    for xi in Base.OneTo(Nx)
    #        ehkslice = ehk2[[xi, xi+Nx], :]
    #        #前向
    #        funcforward = (n) -> bmat_IsingADX(ehkslice, Val(hscfg[bi, xi]), U, 0.0, n[1], n[2])
    #        jacxi = jacobian(funcforward, [ny[xi], nz[xi]])
    #        #这边不同tau时候的操作其实是相同的，应该可以
    #        #[11(hscfg[1,xi],n) 12(hscfg[1,xi],n); 11(hscfg[tau,xi],n) 12(hscfg[tau,xi],n)]
    #        #1. [M11(1) M12(1); ...; M11(tau) M12(tau)] * ehkslice -> (tau, 2*Nx)
    #        #2. [M21(1) M22(1); ...; M21(tau) M22(tau)] * ehkslice -> (tau, 2*Nx)
    #        #reshape(1., tau, 1, 2*Nx) ...
    #        #cat(real1, real2, imag1, imag2; dims=2) -> (tau, 4, 2*Nx)
    #        #diff -> (tau, 4, 2*Nx, 2)
    #        #@time jacxi = @autojac(ehkslice, hscfg, bi, xi, U, ny, nz)
    #        #println("_")
    #        #@time jacxi_ = @autojac2(ehkslice, hscfg, bi, xi, U, ny, nz)
    #        #println(size(jacxi), size(jacxi_))
    #        #@assert all(isapprox.(jacxi, jacxi_))
    #        jacxi = reshape(jacxi, 4, 2*Nx, 2)
    #        pbrpny[bi, [xi, xi+Nx], :, xi] = jacxi[1:2, :, 1]
    #        pbrpnz[bi, [xi, xi+Nx], :, xi] = jacxi[1:2, :, 2]
    #        pbipny[bi, [xi, xi+Nx], :, xi] = jacxi[3:4, :, 1]
    #        pbipnz[bi, [xi, xi+Nx], :, xi] = jacxi[3:4, :, 2]
    #    end
    #    #end
    #    #println(pbrpny[bi, :, :, :], b1)
    #    #println(pbrpnz[bi, :, :, :], b2)
    #    #println(maximum(pbrpny[bi, :, :, :] - b1), maximum(pbrpnz[bi, :, :, :] - b2))
    #    #@assert all(isapprox.(pbrpny[bi, :, :, :], b1))
    #    #@assert all(isapprox.(pbrpnz[bi, :, :, :], b2))
    #    #println(maximum(pbipny[bi, :, :, :] - b3), maximum(pbipnz[bi, :, :, :] - b4))
    #    #@assert all(isapprox.(pbipny[bi, :, :, :], b3))
    #    #@assert all(isapprox.(pbipnz[bi, :, :, :], b4))
    #    #pbrpny[bi, :, :, :] = reshape(brny, 2*Nx, 2*Nx, Nx)
    #    #pbrpnz[bi, :, :, :] = reshape(brnz, 2*Nx, 2*Nx, Nx)
    #    #
    #    #@time jacobian!(result, jfbi, (hscfg_f64[bi, :], ny, nz))
    #    #_, biny, binz = result
    #    #b1 = biny
    #    #@time _, biny, binz = jacobian(fbi, hscfg_f64[bi, :], ny, nz)
    #    #b2 = biny
    #    #println(size(b1), size(b2))
    #    #@assert all(isapprox.(reshape(b1, 2*Nx, 2*Nx, Nx), reshape(b2, 2*Nx, 2*Nx, Nx), atol=1e-13))
    #    #@time begin
    #    #for xi in Base.OneTo(Nx)
    #    #    ehkslice = ehk2[[xi, xi+Nx], :]
    #    #    #前向
    #    #    funcforward = (n) -> bmat_IsingADX(ehkslice, Val(hscfg[bi, xi]), U, 0.0, n[1], n[2])[2]
    #    #    jacxi = jacobian(funcforward, [ny[xi], nz[xi]])
    #    #    pbipny[bi, [xi, xi+Nx], :, xi] = reshape(jacxi[:, 1], 2, 2*Nx)
    #    #    pbipnz[bi, [xi, xi+Nx], :, xi] = reshape(jacxi[:, 2], 2, 2*Nx)
    #    #end
    #    #end
    #    #pbipny[bi, :, :, :] = reshape(biny, 2*Nx, 2*Nx, Nx)
    #    #pbipnz[bi, :, :, :] = reshape(binz, 2*Nx, 2*Nx, Nx)
    #    #
    #end
    #end
    #println("tau")
    #这里内存用的太多了，优化成随bi和xi走的，节约内存
    pbrpnx2 = zeros(Nt, 2*Nx, 2*Nx, Nx)#Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbipnx2 = zeros(Nt, 2*Nx, 2*Nx, Nx)
    pbrpny2 = zeros(Nt, 2*Nx, 2*Nx, Nx)#Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbipny2 = zeros(Nt, 2*Nx, 2*Nx, Nx)#Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbrpnz2 = zeros(Nt, 2*Nx, 2*Nx, Nx)#Array{Float64, 4}(undef, Nt, 2*Nx, 2*Nx, Nx)
    pbipnz2 = zeros(Nt, 2*Nx, 2*Nx, Nx)
    #@time begin
    for xi in Base.OneTo(Nx)
        ehkslice = ehk2[[xi, xi+Nx], :]
        funcforward2 = (n) -> bmat_IsingADTX(ehkslice, hscfg[:, xi], U, n[1], n[2], n[3])
        jacxit = jacobian(funcforward2, [nx[xi], ny[xi], nz[xi]])
        jacxit = reshape(jacxit, 4*Nt, 2*Nx, 3)
        #funcforward3 = (ehk, ths, n) ->  bmat_IsingADTX(ehk, ths, U, 0.0, n[1], n[2])
        #println("Enzyme")
        #@time jacxit2 = cat(autodiff(Forward, funcforward3, BatchDuplicated, Const(ehkslice), Const(hscfg[:, xi]),
        #BatchDuplicated([ny[xi], nz[xi]], ([1.0, 0.0], [0.0, 1.0])))[2]...; 
        #dims=3)
        #println(size(jacxit2))
        #@assert all(isapprox.(jacxit, jacxit2))
        pbrpnx2[:, xi, :, xi]    = jacxit[1:Nt, :, 1]
        pbrpnx2[:, xi+Nx, :, xi] = jacxit[(Nt+1):2Nt, :, 1]
        pbrpny2[:, xi, :, xi]    = jacxit[1:Nt, :, 2]
        pbrpny2[:, xi+Nx, :, xi] = jacxit[(Nt+1):2Nt, :, 2]
        pbrpnz2[:, xi, :, xi]    = jacxit[1:Nt, :, 3]
        pbrpnz2[:, xi+Nx, :, xi] = jacxit[(Nt+1):2Nt, :, 3]
        #
        pbipnx2[:, xi, :, xi]    = jacxit[(2Nt+1):3Nt, :, 1]
        pbipnx2[:, xi+Nx, :, xi] = jacxit[(3Nt+1):4Nt, :, 1]
        pbipny2[:, xi, :, xi]    = jacxit[(2Nt+1):3Nt, :, 2]
        pbipny2[:, xi+Nx, :, xi] = jacxit[(3Nt+1):4Nt, :, 2]
        pbipnz2[:, xi, :, xi]    = jacxit[(2Nt+1):3Nt, :, 3]
        pbipnz2[:, xi+Nx, :, xi] = jacxit[(3Nt+1):4Nt, :, 3]
        #println(jacxit[1, 1:2, :, 1], pbrpny[1, [1, 1+Nx], :, 1])
        #@assert all(isapprox.(jacxit[:, 1:2, :, 1], pbrpny[:, [xi, xi+Nx], :, xi]))
        #@assert all(isapprox.(jacxit[:, 1:2, :, 2], pbrpnz[:, [xi, xi+Nx], :, xi]))
        #@assert all(isapprox.(jacxit[:, 3:4, :, 1], pbipny[:, [xi, xi+Nx], :, xi]))
        #@assert all(isapprox.(jacxit[:, 3:4, :, 2], pbipnz[:, [xi, xi+Nx], :, xi]))
    end
    #end
    return pbrpnx2, pbipnx2, pbrpny2, pbipny2, pbrpnz2, pbipnz2
end
=#

"""
计算每个格点上的jactape，在计算∂B(τ)/∂n_{x,y,z;i}的时候，所有的i之间和τ之间没有什么关联，
每个格点和虚时都单看只有辅助场s_{τ,i}=±1的两种，把他们提前算出来，然后在不同HS场下平均的时候就不要重复计算了。
"""
function pbpn_calc_tapemap(sp::Splitting, nx::Vector{Float64}, ny::Vector{Float64}, nz::Vector{Float64})
    #
    Nx = length(nx)
    tapemap = Matrix{AutoJacTape}(undef, Nx, 2*sp.segl)
    for si in Base.OneTo(sp.segl)
        ehk2 = reshape(sp.ehkt[si, :], sp.sizehk)
        Ui = sp.Uit[si, :]
        for xi in Base.OneTo(Nx)
            ehkslice = ehk2[[xi, xi+Nx], :]
            if abs(Ui[xi]) > 1e-6
                funcforward1 = (n) -> bmat_IsingADX(ehkslice, Val(-1), Ui[xi], n[1], n[2], n[3], sp.dt)
                funcforward2 = (n) -> bmat_IsingADX(ehkslice, Val(1), Ui[xi], n[1], n[2], n[3], sp.dt)
                jacxi1 = jacobian(funcforward1, [nx[xi], ny[xi], nz[xi]])
                jacxi2 = jacobian(funcforward2, [nx[xi], ny[xi], nz[xi]])
            else
                jacxi1 = zeros(Float64, 4*2*Nx*3)
                jacxi2 = zeros(Float64, 4*2*Nx*3)
            end
            tapemap[xi, 2*si-1] = AutoJacTape([nx[xi], ny[xi], nz[xi]], jacxi1; ehkslice=ehkslice, val=-1, U=Ui[xi])
            tapemap[xi, 2*si] = AutoJacTape([nx[xi], ny[xi], nz[xi]], jacxi2; ehkslice=ehkslice, val=1, U=Ui[xi])
        end
    end
    return tapemap
end


"""
计算B矩阵的正向导数
"""
function pbpn_calc_xi(sp, hsxi::Vector{Int}, xidx::Int, nxi::Float64, nyi::Float64, nzi::Float64;
    tapemap=missing)
    #
    Nt = length(hsxi)
    Nx = length(sp.Uit[1, :])
    pbrpnx2 = zeros(Nt, 2, 2*Nx)
    pbipnx2 = zeros(Nt, 2, 2*Nx)
    pbrpny2 = zeros(Nt, 2, 2*Nx)
    pbipny2 = zeros(Nt, 2, 2*Nx)
    pbrpnz2 = zeros(Nt, 2, 2*Nx)
    pbipnz2 = zeros(Nt, 2, 2*Nx)
    #funcforward2 = (n) -> bmat_IsingADTX(ehkslice, hsxi, Ui, n[1], n[2], n[3])
    #jacxit = jacobian(funcforward2, [nxi, nyi, nzi])
    #jacxit = reshape(jacxit, 4*Nt, 2*Nx, 3)
    #
    if !ismissing(tapemap)
        for bi in Base.OneTo(Nt)
            #jac1 = vcat(jacxit[[bi, Nt+bi], :, :], jacxit[[2Nt+bi, 3Nt+bi], :, :])
            segl = sp.vecl[bi]
            tape = hsxi[bi] == 1 ? tapemap[xidx, 2*segl] : tapemap[xidx, 2*segl-1]
            jac1 = reshape(tape.output, 4, 2*Nx, 3)
            #@assert all(isapprox.(tape.input, [nxi, nyi, nzi]))
            #@assert all(isapprox.(jac1, reshape(tape.output, 4, 2*Nx, 3)))
            pbrpnx2[bi, :, :] = jac1[1:2, :, 1]
            pbrpny2[bi, :, :] = jac1[1:2, :, 2]
            pbrpnz2[bi, :, :] = jac1[1:2, :, 3]
            pbipnx2[bi, :, :] = jac1[3:4, :, 1]
            pbipny2[bi, :, :] = jac1[3:4, :, 2]
            pbipnz2[bi, :, :] = jac1[3:4, :, 3]
        end
    else
        @warn "计算jacobian"
        #Ui = sp.Uit[:, xidx]
        #Ui = sum(Ui) / length(Ui)
        for bi in Base.OneTo(Nt)
            ehk2, Ui = unpack_splitting(sp, bi)
            funcforward2 = (n) -> bmat_IsingADX(ehk2[[xidx,xidx+Nx], :],
            Val(hsxi[bi]), Ui[xidx], n[1], n[2], n[3], sp.dt)
            jacxit = jacobian(funcforward2, [nxi, nyi, nzi])
            jac1 = reshape(jacxit, 4, 2*Nx, 3)
            pbrpnx2[bi, :, :] = jac1[1:2, :, 1]
            pbrpny2[bi, :, :] = jac1[1:2, :, 2]
            pbrpnz2[bi, :, :] = jac1[1:2, :, 3]
            pbipnx2[bi, :, :] = jac1[3:4, :, 1]
            pbipny2[bi, :, :] = jac1[3:4, :, 2]
            pbipnz2[bi, :, :] = jac1[3:4, :, 3]
        end
        #funcforward2 = (n) -> bmat_IsingADTX(ehkslice, hsxi, Ui, n[1], n[2], n[3])
        #jacxit = jacobian(funcforward2, [nxi, nyi, nzi])
        #jacxit = reshape(jacxit, 4*Nt, 2*Nx, 3)
        #pbrpnx2[:, 1, :] = jacxit[1:Nt, :, 1]
        #pbrpnx2[:, 2, :] = jacxit[(Nt+1):2Nt, :, 1]
        #pbrpny2[:, 1, :] = jacxit[1:Nt, :, 2]
        #pbrpny2[:, 2, :] = jacxit[(Nt+1):2Nt, :, 2]
        #pbrpnz2[:, 1, :] = jacxit[1:Nt, :, 3]
        #pbrpnz2[:, 2, :] = jacxit[(Nt+1):2Nt, :, 3]
        ##
        #pbipnx2[:, 1, :] = jacxit[(2Nt+1):3Nt, :, 1]
        #pbipnx2[:, 2, :] = jacxit[(3Nt+1):4Nt, :, 1]
        #pbipny2[:, 1, :] = jacxit[(2Nt+1):3Nt, :, 2]
        #pbipny2[:, 2, :] = jacxit[(3Nt+1):4Nt, :, 2]
        #pbipnz2[:, 1, :] = jacxit[(2Nt+1):3Nt, :, 3]
        #pbipnz2[:, 2, :] = jacxit[(3Nt+1):4Nt, :, 3]
    end
    return pbrpnx2, pbipnx2, pbrpny2, pbipny2, pbrpnz2, pbipnz2
end


"""
计算B矩阵的正向导数
"""
function pbpn_calc(sp, hscfg::Matrix{Int}, nx::Vector{Float64}, ny::Vector{Float64}, nz::Vector{Float64};
    tapemap=missing)
    Nt, Nx = size(hscfg)
    #
    pbrpnx2 = zeros(Nt, 2*Nx, 2*Nx, Nx)
    pbipnx2 = zeros(Nt, 2*Nx, 2*Nx, Nx)
    pbrpny2 = zeros(Nt, 2*Nx, 2*Nx, Nx)
    pbipny2 = zeros(Nt, 2*Nx, 2*Nx, Nx)
    pbrpnz2 = zeros(Nt, 2*Nx, 2*Nx, Nx)
    pbipnz2 = zeros(Nt, 2*Nx, 2*Nx, Nx)
    for xi in Base.OneTo(Nx)
        rpnx, ipnx, rpny, ipny, rpnz, ipnz = pbpn_calc_xi(
            sp, hscfg[:, xi], xi, nx[xi], ny[xi], nz[xi];
            tapemap=tapemap
        )
        pbrpnx2[:, [xi, xi+Nx], :, xi] = rpnx
        pbipnx2[:, [xi, xi+Nx], :, xi] = ipnx
        pbrpny2[:, [xi, xi+Nx], :, xi] = rpny
        pbipny2[:, [xi, xi+Nx], :, xi] = ipny
        pbrpnz2[:, [xi, xi+Nx], :, xi] = rpnz
        pbipnz2[:, [xi, xi+Nx], :, xi] = ipnz
    end
    return pbrpnx2, pbipnx2, pbrpny2, pbipny2, pbrpnz2, pbipnz2
end


"""
计算∂L/∂n
"""
function meas_grad(ss::ScrollSVD{T}, allbmats, sp, allcfgs, nx, ny, nz;
    tapemap=missing) where T
    Nt = length(allbmats)
    Nx = length(nx)
    # 计算 \bar{G} = (1 + B^-1)^-1
    # 计算 \bar{B_l}_{xy} = -(GB)^T
    bbar = bbar_calc(ss, allbmats)
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{theta_i} = \sum \bar{B_l}_xy  ∂B / ∂theta_i
    nxbar2 = zeros(Float64, Nx)
    nybar2 = zeros(Float64, Nx)
    nzbar2 = zeros(Float64, Nx)
    for bi in Base.OneTo(Nt)
        ehk2, Ui = unpack_splitting(sp, bi)
        for xi in Base.OneTo(Nx)
            ehkslice = ehk2[[xi, xi+Nx], :]
            rpnx, ipnx, rpny, ipny, rpnz, ipnz = pbpn_calc_xi(
                sp, allcfgs[:, xi], xi, nx[xi], ny[xi], nz[xi];
                tapemap=tapemap
            )
            #注意虚部附加的负号
            #xi只会和xi，xi+Nx能向前传递
            nzbar2[xi] += sum(rpnz[bi, 1, :] .* real(bbar[bi, xi, :]))
            nzbar2[xi] += sum(rpnz[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            nzbar2[xi] -= sum(ipnz[bi, 1, :] .* imag(bbar[bi, xi, :]))
            nzbar2[xi] -= sum(ipnz[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
            nybar2[xi] += sum(rpny[bi, 1, :] .* real(bbar[bi, xi, :]))
            nybar2[xi] += sum(rpny[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            nybar2[xi] -= sum(ipny[bi, 1, :] .* imag(bbar[bi, xi, :]))
            nybar2[xi] -= sum(ipny[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
            nxbar2[xi] += sum(rpnx[bi, 1, :] .* real(bbar[bi, xi, :]))
            nxbar2[xi] += sum(rpnx[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            nxbar2[xi] -= sum(ipnx[bi, 1, :] .* imag(bbar[bi, xi, :]))
            nxbar2[xi] -= sum(ipnx[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
        end
    end
    return nxbar2, nybar2, nzbar2
end


#=
"""
计算∂L/∂n
"""
function meas_grad_(ss::ScrollSVD{T}, allbmats, hk, allcfgs, Ui, nx, ny, nz; Z2=true) where T
    Nt = length(allbmats)
    Nx = length(Ui)
    # 计算 \bar{G} = (1 + B^-1)^-1
    # 计算 \bar{B_l}_{xy} = -(GB)^T
    bbar = bbar_calc(ss, allbmats)
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    Uval = sum(Ui) / length(Ui)
    pbrpnx, pbipnx, pbrpny, pbipny, pbrpnz, pbipnz = pbpn_calc(hk, allcfgs, Uval, nx, ny, nz; Z2=Z2)
    dt = 0.1
    ehk = exp(-dt*hk)
    if Z2
        ehk2 = kron([1 0; 0 1], ehk)
    else
        ehk2 = ehk
    end
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{theta_i} = \sum \bar{B_l}_xy  ∂B / ∂theta_i
    nxbar2 = zeros(Float64, Nx)
    nybar2 = zeros(Float64, Nx)
    nzbar2 = zeros(Float64, Nx)
    for xi in Base.OneTo(Nx)
        ehkslice = ehk2[[xi, xi+Nx], :]
        rpnx, ipnx, rpny, ipny, rpnz, ipnz = pbpn_calc_xi(
            ehkslice, allcfgs[:, xi], Ui[xi], nx[xi], ny[xi], nz[xi]
        )
        #@assert all(isapprox.(pbrpnx[:, [xi, xi+Nx], :, xi], rpnx))
        #@assert all(isapprox.(pbipnx[:, [xi, xi+Nx], :, xi], ipnx))
        #@assert all(isapprox.(pbrpny[:, [xi, xi+Nx], :, xi], rpny))
        #@assert all(isapprox.(pbipny[:, [xi, xi+Nx], :, xi], ipny))
        #@assert all(isapprox.(pbrpnz[:, [xi, xi+Nx], :, xi], rpnz))
        #@assert all(isapprox.(pbipnz[:, [xi, xi+Nx], :, xi], ipnz))
        #xi只会和xi，xi+Nx能向前传递
        nzbar2[xi] += sum(rpnz[:, 1, :] .* real(bbar[:, xi, :]))
        nzbar2[xi] += sum(rpnz[:, 2, :] .* real(bbar[:, xi+Nx, :]))
        nzbar2[xi] -= sum(ipnz[:, 1, :] .* imag(bbar[:, xi, :]))
        nzbar2[xi] -= sum(ipnz[:, 2, :] .* imag(bbar[:, xi+Nx, :]))
        nybar2[xi] += sum(rpny[:, 1, :] .* real(bbar[:, xi, :]))
        nybar2[xi] += sum(rpny[:, 2, :] .* real(bbar[:, xi+Nx, :]))
        nybar2[xi] -= sum(ipny[:, 1, :] .* imag(bbar[:, xi, :]))
        nybar2[xi] -= sum(ipny[:, 2, :] .* imag(bbar[:, xi+Nx, :]))
        nxbar2[xi] += sum(rpnx[:, 1, :] .* real(bbar[:, xi, :]))
        nxbar2[xi] += sum(rpnx[:, 2, :] .* real(bbar[:, xi+Nx, :]))
        nxbar2[xi] -= sum(ipnx[:, 1, :] .* imag(bbar[:, xi, :]))
        nxbar2[xi] -= sum(ipnx[:, 2, :] .* imag(bbar[:, xi+Nx, :]))
    end
    #
    nxbar = zeros(Float64, Nx)
    nybar = zeros(Float64, Nx)
    nzbar = zeros(Float64, Nx)
    #注意虚部附加的负号
    #TODO：这里先算出来整个的pbrpnz实在太消耗内存了，优化成计算某个时刻和位置的
    for bi in Base.OneTo(Nt); for xi in Base.OneTo(2*Nx); for yi in Base.OneTo(2*Nx)
        nzbar += pbrpnz[bi, xi, yi, :] * real(bbar[bi, xi, yi])
        nzbar -= pbipnz[bi, xi, yi, :] * imag(bbar[bi, xi, yi])
        nybar += pbrpny[bi, xi, yi, :] * real(bbar[bi, xi, yi])
        nybar -= pbipny[bi, xi, yi, :] * imag(bbar[bi, xi, yi])
        nxbar += pbrpnx[bi, xi, yi, :] * real(bbar[bi, xi, yi])
        nxbar -= pbipnx[bi, xi, yi, :] * imag(bbar[bi, xi, yi])
    end; end; end
    @assert all(isapprox.(nxbar, nxbar2))
    @assert all(isapprox.(nybar, nybar2))
    @assert all(isapprox.(nzbar, nzbar2))
    return nxbar2, nybar2, nzbar2
end
=#



"""
计算每个格点上的∂B(τ)/∂a
"""
function pbpa_calc_tapemap(sp::Splitting, rax::Vector{Float64}, iax::Vector{Float64}, rbx::Vector{Float64}, ibx::Vector{Float64})
    #
    Nx = length(rax)
    tapemap = Matrix{AutoJacTape}(undef, Nx, 2*sp.segl)
    for si in Base.OneTo(sp.segl)
        ehk2 = reshape(sp.ehkt[si, :], sp.sizehk)
        Ui = sp.Uit[si, :]
        for xi in Base.OneTo(Nx)
            ehkslice = ehk2[[xi, xi+Nx], :]
            if abs(Ui[xi]) > 1e-6
                funcforward1 = (n) -> bmat_Gauge1ADX(ehkslice, Val(-1), Ui[xi], n[1], n[2], n[3], n[4], sp.dt)
                funcforward2 = (n) -> bmat_Gauge1ADX(ehkslice, Val(1), Ui[xi], n[1], n[2], n[3], n[4], sp.dt)
                jacxi1 = jacobian(funcforward1, [rax[xi], iax[xi], rbx[xi], ibx[xi]])
                jacxi2 = jacobian(funcforward2, [rax[xi], iax[xi], rbx[xi], ibx[xi]])
            else
                jacxi1 = zeros(Float64, 4*2*Nx*4)
                jacxi2 = zeros(Float64, 4*2*Nx*4)
            end
            tapemap[xi, 2*si-1] = AutoJacTape([rax[xi], iax[xi], rbx[xi], ibx[xi]], jacxi1; val=-1, U=Ui[xi])
            tapemap[xi, 2*si] = AutoJacTape([rax[xi], iax[xi], rbx[xi], ibx[xi]], jacxi2; val=1, U=Ui[xi])
        end
    end
    return tapemap
end


"""
计算B矩阵的正向导数
"""
function pbpa_calc_xi(sp, hsxi::Vector{Int}, xidx::Int, ra::Float64, ia::Float64, rb::Float64, ib::Float64;
    tapemap=missing)
    #
    Nt = length(hsxi)
    Nx = length(sp.Uit[1, :])
    pbrpra2 = zeros(Nt, 2, 2*Nx)
    pbipra2 = zeros(Nt, 2, 2*Nx)
    pbrpia2 = zeros(Nt, 2, 2*Nx)
    pbipia2 = zeros(Nt, 2, 2*Nx)
    pbrprb2 = zeros(Nt, 2, 2*Nx)
    pbiprb2 = zeros(Nt, 2, 2*Nx)
    pbrpib2 = zeros(Nt, 2, 2*Nx)
    pbipib2 = zeros(Nt, 2, 2*Nx)
    #
    if !ismissing(tapemap)
        for bi in Base.OneTo(Nt)
            #jac1 = vcat(jacxit[[bi, Nt+bi], :, :], jacxit[[2Nt+bi, 3Nt+bi], :, :])
            segl = sp.vecl[bi]
            tape = hsxi[bi] == 1 ? tapemap[xidx, 2*segl] : tapemap[xidx, 2*segl-1]
            jac1 = reshape(tape.output, 4, 2*Nx, 4)
            #B的实数部分在1:2
            pbrpra2[bi, :, :] = jac1[1:2, :, 1]
            pbrpia2[bi, :, :] = jac1[1:2, :, 2]
            pbrprb2[bi, :, :] = jac1[1:2, :, 3]
            pbrpib2[bi, :, :] = jac1[1:2, :, 4]
            #B的虚数部分在3:4
            pbipra2[bi, :, :] = jac1[3:4, :, 1]
            pbipia2[bi, :, :] = jac1[3:4, :, 2]
            pbiprb2[bi, :, :] = jac1[3:4, :, 3]
            pbipib2[bi, :, :] = jac1[3:4, :, 4]
        end
    else
        @warn "计算jacobian"
        #Ui = sp.Uit[:, xidx]
        #Ui = sum(Ui) / length(Ui)
        for bi in Base.OneTo(Nt)
            ehk2, Ui = unpack_splitting(sp, bi)
            funcforward2 = (n) -> bmat_Gauge1ADX(ehk2[[xidx,xidx+Nx], :],
            Val(hsxi[bi]), Ui[xidx], n[1], n[2], n[3], n[4], sp.dt)
            jacxit = jacobian(funcforward2, [ra, ia, rb, ib])
            jac1 = reshape(jacxit, 4, 2*Nx, 4)
            pbrpra2[bi, :, :] = jac1[1:2, :, 1]
            pbrpia2[bi, :, :] = jac1[1:2, :, 2]
            pbrprb2[bi, :, :] = jac1[1:2, :, 3]
            pbrpib2[bi, :, :] = jac1[1:2, :, 4]
            pbipra2[bi, :, :] = jac1[3:4, :, 1]
            pbipia2[bi, :, :] = jac1[3:4, :, 2]
            pbiprb2[bi, :, :] = jac1[3:4, :, 3]
            pbipib2[bi, :, :] = jac1[3:4, :, 4]
        end
    end
    return pbrpra2, pbipra2, pbrpia2, pbipia2, pbrprb2, pbiprb2, pbrpib2, pbipib2
end


"""
计算∂L/∂a
"""
function meas_grada(ss::ScrollSVD{T}, allbmats, sp, allcfgs, rax, iax, rbx, ibx;
    tapemap=missing, ignorew=false) where T
    Nt, Nx = size(allcfgs)
    # 计算 \bar{G} = (1 + B^-1)^-1
    # 计算 \bar{B_l}_{xy} = -(GB)^T
    bbar = bbar_calc(ss, allbmats)
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{theta_i} = \sum \bar{B_l}_xy  ∂B / ∂theta_i
    rabar2 = zeros(Float64, Nx)
    iabar2 = zeros(Float64, Nx)
    rbbar2 = zeros(Float64, Nx)
    ibbar2 = zeros(Float64, Nx)
    pval = @. rbx + im*ibx
    for bi in Base.OneTo(Nt)
        #ehk2, Ui = unpack_splitting(sp, bi)
        for xi in Base.OneTo(Nx)
            #ehkslice = ehk2[[xi, xi+Nx], :]
            rpra, ipra, rpia, ipia, rprb, iprb, rpib, ipib = pbpa_calc_xi(
                sp, allcfgs[:, xi], xi, rax[xi], iax[xi], rbx[xi], ibx[xi];
                tapemap=tapemap
            )
            #注意虚部附加的负号
            #xi只会和xi，xi+Nx能向前传递
            rabar2[xi] += sum(rpra[bi, 1, :] .* real(bbar[bi, xi, :]))
            rabar2[xi] += sum(rpra[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            rabar2[xi] -= sum(ipra[bi, 1, :] .* imag(bbar[bi, xi, :]))
            rabar2[xi] -= sum(ipra[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
            #
            iabar2[xi] += sum(rpia[bi, 1, :] .* real(bbar[bi, xi, :]))
            iabar2[xi] += sum(rpia[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            iabar2[xi] -= sum(ipia[bi, 1, :] .* imag(bbar[bi, xi, :]))
            iabar2[xi] -= sum(ipia[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
            #
            rbbar2[xi] += sum(rprb[bi, 1, :] .* real(bbar[bi, xi, :]))
            rbbar2[xi] += sum(rprb[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            rbbar2[xi] -= sum(iprb[bi, 1, :] .* imag(bbar[bi, xi, :]))
            rbbar2[xi] -= sum(iprb[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
            #
            ibbar2[xi] += sum(rpib[bi, 1, :] .* real(bbar[bi, xi, :]))
            ibbar2[xi] += sum(rpib[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            ibbar2[xi] -= sum(ipib[bi, 1, :] .* imag(bbar[bi, xi, :]))
            ibbar2[xi] -= sum(ipib[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
            #
            if ignorew
                continue
            end
            #
            #if allcfgs[bi, xi] == -1
            #wgt = rb + i*ib
            #-Re ∂lnwgt/∂rb = -Re{1 / (rb + i*ib)} = -rb / (rb^2 + ib^2)
            #-Re ∂lnwgt/∂ib = -Re{i / (rb + i*ib)} = -ib / (rb^2 + ib^2)
            #else
            #wgt = 1 -rb -i*ib
            #-Re ∂lnwgt/∂rb = -Re{-1 / (1 -rb -i*ib)} = (1-rb) / ((1-rb)^2 + ib^2)
            #-Re ∂lnwgt/∂ri = -Re{-i / (1 -rb -i*ib)} = -ib / ((1-rb)^2 + ib^2)
            #注意dlnS = dlnZ - dlnw，这个负号在B矩阵时是包含在bbar中的
            #现在包含在导数的结果里
            if allcfgs[bi, xi] == -1
                rbbar2[xi] += -rbx[xi] / (rbx[xi]^2 + ibx[xi]^2)
                ibbar2[xi] += -ibx[xi] / (rbx[xi]^2 + ibx[xi]^2)
            else
                rbbar2[xi] += (1-rbx[xi]) / ((1-rbx[xi])^2 + ibx[xi]^2)
                ibbar2[xi] += -ibx[xi] / ((1-rbx[xi])^2 + ibx[xi]^2)
            end
            #    Uibar2[xi] -= sqrtma*ηdict[allcfgs[bi, xi]]
        end
    end
    return rabar2, iabar2, rbbar2, ibbar2
end


"""
计算每个格点上的jactape，在计算∂B(τ)/∂U_{i}的时候，所有的i之间和τ之间没有什么关联，
每个格点和虚时都单看只有辅助场s_{τ,i}=±1的两种，把他们提前算出来，然后在不同HS场下平均的时候就不要重复计算了。
"""
function pbpU_calc_tapemap(sp::Splitting)
    #
    Nx = size(sp.Uit)[2]
    tapemap = Matrix{AutoJacTape}(undef, Nx, 2*sp.segl)
    for si in Base.OneTo(sp.segl)
        ehk2 = reshape(sp.ehkt[si, :], sp.sizehk)
        Ui = sp.Uit[si, :]
        for xi in Base.OneTo(Nx)
            ehkslice = ehk2[[xi, xi+Nx], :]
            if abs(Ui[xi]) > 1e-6
                funcforward1 = (U) -> bmat_IsingADX(ehkslice, Val(-1), U, 0.0, 0.0, 1.0, sp.dt)
                funcforward2 = (U) -> bmat_IsingADX(ehkslice, Val(1), U, 0.0, 0.0, 1.0, sp.dt)
                jacxi1 = derivative(funcforward1, Ui[xi])
                jacxi2 = derivative(funcforward2, Ui[xi])
            else
                jacxi1 = zeros(Float64, 4*2*Nx)
                jacxi2 = zeros(Float64, 4*2*Nx)
            end
            tapemap[xi, 2*si-1] = AutoJacTape([Ui[xi]], jacxi1; val=-1, U=Ui[xi])
            tapemap[xi, 2*si] = AutoJacTape([Ui[xi]], jacxi2; val=1, U=Ui[xi])
        end
    end
    return tapemap
end


"""
计算B矩阵的正向导数
"""
function pbpU_calc_xi(sp::Splitting, hsxi::Vector{Int}, xidx::Int; tapemap=missing)
    #
    Nt = length(hsxi)
    Nx = length(sp.Uit[1, :])
    pbrpU2 = zeros(Nt, 2, 2*Nx)
    pbipU2 = zeros(Nt, 2, 2*Nx)
    #funcforward2 = (n) -> bmat_IsingADTX(ehkslice, hsxi, Ui, n[1], n[2], n[3])
    #jacxit = jacobian(funcforward2, [nxi, nyi, nzi])
    #jacxit = reshape(jacxit, 4*Nt, 2*Nx, 3)
    #
    if !ismissing(tapemap)
        for bi in Base.OneTo(Nt)
            #jac1 = vcat(jacxit[[bi, Nt+bi], :, :], jacxit[[2Nt+bi, 3Nt+bi], :, :])
            segl = sp.vecl[bi]
            tape = hsxi[bi] == 1 ? tapemap[xidx, 2*segl] : tapemap[xidx, 2*segl-1]
            jac1 = reshape(tape.output, 4, 2*Nx)
            #@assert all(isapprox.(tape.input, [nxi, nyi, nzi]))
            #@assert all(isapprox.(jac1, reshape(tape.output, 4, 2*Nx, 3)))
            pbrpU2[bi, :, :] = jac1[1:2, :]
            pbipU2[bi, :, :] = jac1[3:4, :]
        end
    else
        @warn "计算jacobian"
        #Ui = sp.Uit[:, xidx]
        #Ui = sum(Ui) / length(Ui)
        for bi in Base.OneTo(Nt)
            ehk2, Ui = unpack_splitting(sp, bi)
            #因为bmat_IsingADX中不接受vector的输入，jacibian不会重新定义exp（只有exp.）
            #这里要用derivative
            funcforward2 = (U) -> bmat_IsingADX(ehk2[[xidx,xidx+Nx], :],
            Val(hsxi[bi]), U, 0.0, 0.0, 1.0, sp.dt)
            #
            jacxit = derivative(funcforward2, Ui[xidx])
            jac1 = reshape(jacxit, 4, 2*Nx)
            pbrpU2[bi, :, :] = jac1[1:2, :]
            pbipU2[bi, :, :] = jac1[3:4, :]
        end
    end
    return pbrpU2, pbipU2
end


"""
计算lnS对U的导数
"""
function meas_gradU(ss::ScrollSVD{T}, allbmats, sp, allcfgs; tapemap=missing) where T
    Nt = length(allbmats)
    Nx = size(sp.Uit)[2]
    # 计算 \bar{G} = (1 + B^-1)^-1
    # 计算 \bar{B_l}_{xy} = -(GB)^T
    bbar = bbar_calc(ss, allbmats)
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{theta_i} = \sum \bar{B_l}_xy  ∂B / ∂theta_i
    Uibar2 = zeros(Float64, Nx)
    for xi in Base.OneTo(Nx)
        rpU, ipU = pbpU_calc_xi(
            sp, allcfgs[:, xi], xi;tapemap=tapemap
        )
        for bi in Base.OneTo(Nt)
            #注意虚部附加的负号
            #xi只会和xi，xi+Nx能向前传递
            Uibar2[xi] += sum(rpU[bi, 1, :] .* real(bbar[bi, xi, :]))
            Uibar2[xi] += sum(rpU[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            Uibar2[xi] -= sum(ipU[bi, 1, :] .* imag(bbar[bi, xi, :]))
            Uibar2[xi] -= sum(ipU[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
        end
    end
    return Uibar2
end


"""
计算每个格点上的jactape，在计算∂B(τ)/∂U_{i}的时候，所有的i之间和τ之间没有什么关联，
每个格点和虚时都单看只有辅助场s_{τ,i}=±1的两种，把他们提前算出来，然后在不同HS场下平均的时候就不要重复计算了。
"""
function pbpU_calc_tapemap_Quad1(sp::Splitting, ϕs::Vector{Float64})
    Nx = size(sp.Uit)[2]
    tapemap = Matrix{AutoJacTape}(undef, Nx, 4*sp.segl)
    for si in Base.OneTo(sp.segl)
        ehk2 = reshape(sp.ehkt[si, :], sp.sizehk)
        Ui = sp.Uit[si, :]
        for xi in Base.OneTo(Nx)
            ehkslice = ehk2[[xi, xi+Nx], :]
            funcforward1 = (u) -> bmat_Quad1ADX(ehkslice, Val(-2), u, ϕs[xi], sp.dt)
            funcforward2 = (u) -> bmat_Quad1ADX(ehkslice, Val(-1), u, ϕs[xi], sp.dt)
            funcforward3 = (u) -> bmat_Quad1ADX(ehkslice, Val(1), u, ϕs[xi], sp.dt)
            funcforward4 = (u) -> bmat_Quad1ADX(ehkslice, Val(2), u, ϕs[xi], sp.dt)
            jacxi1 = derivative(funcforward1, Ui[xi])
            jacxi2 = derivative(funcforward2, Ui[xi])
            jacxi3 = derivative(funcforward3, Ui[xi])
            jacxi4 = derivative(funcforward4, Ui[xi])
            tapemap[xi, 4*si-3] = AutoJacTape([Ui[xi]], jacxi1; val=-2, Ui=Ui[xi])
            tapemap[xi, 4*si-2] = AutoJacTape([Ui[xi]], jacxi2; val=-1, Ui=Ui[xi])
            tapemap[xi, 4*si-1] = AutoJacTape([Ui[xi]], jacxi3; val=1, Ui=Ui[xi])
            tapemap[xi, 4*si] = AutoJacTape([Ui[xi]], jacxi4; val=2, Ui=Ui[xi])
        end
    end
    return tapemap
end


"""
计算B矩阵的正向导数
"""
function pbpU_calc_xi_Quad1(sp::Splitting, hsxi::Vector{Int}, xidx::Int, ϕs::Vector{Float64};
    tapemap=missing)
    #
    Nt = length(hsxi)
    Nx = length(sp.Uit[1, :])
    pbrpU2 = zeros(Nt, 2, 2*Nx)
    pbipU2 = zeros(Nt, 2, 2*Nx)
    #funcforward2 = (n) -> bmat_IsingADTX(ehkslice, hsxi, Ui, n[1], n[2], n[3])
    #jacxit = jacobian(funcforward2, [nxi, nyi, nzi])
    #jacxit = reshape(jacxit, 4*Nt, 2*Nx, 3)
    #
    c2t = Dict(-2 => -3, -1 => -2, 1 => -1, 2 => 0)
    if !ismissing(tapemap)
        for bi in Base.OneTo(Nt)
            #jac1 = vcat(jacxit[[bi, Nt+bi], :, :], jacxit[[2Nt+bi, 3Nt+bi], :, :])
            segl = sp.vecl[bi]
            tape = tapemap[xidx, 4*segl+c2t[hsxi[bi]]]
            jac1 = reshape(tape.output, 4, 2*Nx)
            #@assert all(isapprox.(tape.input, [nxi, nyi, nzi]))
            #@assert all(isapprox.(jac1, reshape(tape.output, 4, 2*Nx, 3)))
            pbrpU2[bi, :, :] = jac1[1:2, :]
            pbipU2[bi, :, :] = jac1[3:4, :]
        end
    else
        @warn "计算jacobian"
        error("not implement")
    end
    return pbrpU2, pbipU2
end


"""
计算lnS对U的导数
"""
function meas_gradU_Quad1(ss::ScrollSVD{T}, allbmats, sp, allcfgs, ϕs;
    tapemap=missing, ignorew=false) where T
    #
    ηdict = Dict(
        -2 => -3.301360247771569, 
        -1 => -1.049295246550581,
        1 => 1.049295246550581,
        2 => 3.301360247771569
    )
    #
    Nt = length(allbmats)
    Nx = size(sp.Uit)[2]
    # 计算 \bar{G} = (1 + B^-1)^-1
    # 计算 \bar{B_l}_{xy} = -(GB)^T
    bbar = bbar_calc(ss, allbmats)
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{theta_i} = \sum \bar{B_l}_xy  ∂B / ∂theta_i
    Uibar2 = zeros(Float64, Nx)
    for xi in Base.OneTo(Nx)
        rpU, ipU = pbpU_calc_xi_Quad1(
            sp, allcfgs[:, xi], xi, ϕs;tapemap=tapemap
        )
        for bi in Base.OneTo(Nt)
            #注意虚部附加的负号
            #xi只会和xi，xi+Nx能向前传递
            Uibar2[xi] += sum(rpU[bi, 1, :] .* real(bbar[bi, xi, :]))
            Uibar2[xi] += sum(rpU[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            Uibar2[xi] -= sum(ipU[bi, 1, :] .* imag(bbar[bi, xi, :]))
            Uibar2[xi] -= sum(ipU[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
        end
        #
        if ignorew
            continue
        end
        #
        for bi in Base.OneTo(Nt)
            #wgt = γ*exp(sqrtma*η*ϕ)
            #注意dlnS = dlnZ - dlnw，这个负号在B矩阵时是包含在bbar中的
            #现在包含在导数的结果里
            #-∂lnwgt/∂U = -(∂lnγ + ∂sqrtma*η*ϕ)/∂U = -η*ϕ*sqrt(dt)/sqrt(8*U)
            Ui = sp.Uit[sp.vecl[bi], xi]
            Uibar2[xi] -= ηdict[allcfgs[bi, xi]]*ϕs[xi]*sqrt(sp.dt)/sqrt(8*Ui)
        end
    end
    return Uibar2
end


"""
计算每个格点上的jactape，在计算∂B(τ)/∂U_{i}的时候，所有的i之间和τ之间没有什么关联，
每个格点和虚时都单看只有辅助场s_{τ,i}=±1的两种，把他们提前算出来，然后在不同HS场下平均的时候就不要重复计算了。
"""
function pbpUabnm_calc_tapemap(sp::Splitting, nx::Vector{Float64}, ny::Vector{Float64}, nz::Vector{Float64})
    #
    Nx = size(sp.Uit)[2]
    tapemap = Matrix{AutoJacTape}(undef, Nx, 2*sp.segl)
    for si in Base.OneTo(sp.segl)
        ehk2 = reshape(sp.ehkt[si, :], sp.sizehk)
        Ui = sp.Uit[si, :]
        for xi in Base.OneTo(Nx)
            ehkslice = ehk2[[xi, xi+Nx], :]
            if abs(Ui[xi]) > 1e-6
                funcforward1 = (U) -> bmat_IsingADX(ehkslice, Val(-1), U, nx[xi], ny[xi], nz[xi], sp.dt)
                funcforward2 = (U) -> bmat_IsingADX(ehkslice, Val(1), U, nx[xi], ny[xi], nz[xi], sp.dt)
                jacxi1 = derivative(funcforward1, Ui[xi])
                jacxi2 = derivative(funcforward2, Ui[xi])
            else
                jacxi1 = zeros(Float64, 4*2*Nx)
                jacxi2 = zeros(Float64, 4*2*Nx)
            end
            tapemap[xi, 2*si-1] = AutoJacTape([Ui[xi]], jacxi1; val=-1, U=Ui[xi])
            tapemap[xi, 2*si] = AutoJacTape([Ui[xi]], jacxi2; val=1, U=Ui[xi])
        end
    end
    return tapemap
end


"""
计算B矩阵的正向导数
"""
function pbpUabnm_calc_xi(sp::Splitting, hsxi::Vector{Int}, xidx::Int,
    nx::Vector{Float64}, ny::Vector{Float64}, nz::Vector{Float64}
    ; tapemap=missing)
    #
    Nt = length(hsxi)
    Nx = length(sp.Uit[1, :])
    pbrpU2 = zeros(Nt, 2, 2*Nx)
    pbipU2 = zeros(Nt, 2, 2*Nx)
    #funcforward2 = (n) -> bmat_IsingADTX(ehkslice, hsxi, Ui, n[1], n[2], n[3])
    #jacxit = jacobian(funcforward2, [nxi, nyi, nzi])
    #jacxit = reshape(jacxit, 4*Nt, 2*Nx, 3)
    #
    if !ismissing(tapemap)
        for bi in Base.OneTo(Nt)
            #jac1 = vcat(jacxit[[bi, Nt+bi], :, :], jacxit[[2Nt+bi, 3Nt+bi], :, :])
            segl = sp.vecl[bi]
            tape = hsxi[bi] == 1 ? tapemap[xidx, 2*segl] : tapemap[xidx, 2*segl-1]
            jac1 = reshape(tape.output, 4, 2*Nx)
            #@assert all(isapprox.(tape.input, [nxi, nyi, nzi]))
            #@assert all(isapprox.(jac1, reshape(tape.output, 4, 2*Nx, 3)))
            pbrpU2[bi, :, :] = jac1[1:2, :]
            pbipU2[bi, :, :] = jac1[3:4, :]
        end
    else
        @warn "计算jacobian"
        #Ui = sp.Uit[:, xidx]
        #Ui = sum(Ui) / length(Ui)
        for bi in Base.OneTo(Nt)
            ehk2, Ui = unpack_splitting(sp, bi)
            #因为bmat_IsingADX中不接受vector的输入，jacibian不会重新定义exp（只有exp.）
            #这里要用derivative
            funcforward2 = (U) -> bmat_IsingADX(ehk2[[xidx,xidx+Nx], :],
            Val(hsxi[bi]), U, nx[xidx], ny[xidx], nz[xidx], sp.dt)
            #
            jacxit = derivative(funcforward2, Ui[xidx])
            jac1 = reshape(jacxit, 4, 2*Nx)
            pbrpU2[bi, :, :] = jac1[1:2, :]
            pbipU2[bi, :, :] = jac1[3:4, :]
        end
    end
    return pbrpU2, pbipU2
end


"""
计算lnS对U的导数
"""
function meas_gradUabnm(ss::ScrollSVD{T}, allbmats, sp, allcfgs, nx, ny, nz; tapemap=missing) where T
    Nt = length(allbmats)
    Nx = size(sp.Uit)[2]
    # 计算 \bar{G} = (1 + B^-1)^-1
    # 计算 \bar{B_l}_{xy} = -(GB)^T
    bbar = bbar_calc(ss, allbmats)
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{n{xyz}_i} = \sum \bar{B_l}_xy  ∂B / ∂U^{xyz}_i 
    # 计算 \bar{theta_i} = \sum \bar{B_l}_xy  ∂B / ∂theta_i
    Uibar2 = zeros(Float64, Nx)
    for xi in Base.OneTo(Nx)
        rpU, ipU = pbpUabnm_calc_xi(
            sp, allcfgs[:, xi], xi, nx, ny, nz;tapemap=tapemap
        )
        for bi in Base.OneTo(Nt)
            #注意虚部附加的负号
            #xi只会和xi，xi+Nx能向前传递
            Uibar2[xi] += sum(rpU[bi, 1, :] .* real(bbar[bi, xi, :]))
            Uibar2[xi] += sum(rpU[bi, 2, :] .* real(bbar[bi, xi+Nx, :]))
            Uibar2[xi] -= sum(ipU[bi, 1, :] .* imag(bbar[bi, xi, :]))
            Uibar2[xi] -= sum(ipU[bi, 2, :] .* imag(bbar[bi, xi+Nx, :]))
        end
    end
    return Uibar2
end
